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We calculate the J/tp mass, leptonic width and radiative decay rate to "/n c from lattice QCD 
including u, d and s quarks in the sea for the first time. We use the Highly Improved Staggered 
Quark formalism and nonperturbatively normalised vector currents for the leptonic and radiative 
decay rates. Our results are: M J/i: - M Vc = 116.5(3.2)MeV; T{J/xj) -> e+e") = 5.48(16)keV; 
r(J/i/> — > 7?7 C ) = 2.49(19)keV. The first two are in good agreement with experiment, with T(J/i/j — >■ 
e + e~) providing a test of a decay matrix element in QCD, independent of CKM uncertainties, to 
2%. At the same time results for the time moments of the correlation function can be compared to 
values from the charm contribution to <r(e + e - — > hadrons), giving a 1.5% test of QCD. Our results 
show that an improved experimental error would enable a similarly strong test from F(J/tp — > yr/c)- 



I. INTRODUCTION 

Precision tests of lattice QCD against experiment are 
critical to provide benchmarks against which to calibrate 
the reliability of predictions from lattice QCD [T] . Most 
tests to date have relied on the spectrum of gold-plated 
hadron masses - for example, the mass of the D s me- 
son can be calculated in lattice QCD with an error of 
3 MeV (having fixed the masses of the c and s masses 
from other mesons) and the result agrees with experi- 
ment 2,3. Here we give another such test by determin- 
ing the mass of the J/ip to a precision of 3 MeV. 

Tests of decay matrix elements are harder to do very 
accurately. We need precision tests of these because it 
is the predictions of decay matrix elements from lattice 
QCD that enable, for example, progress with the flavor 
physics programme [4] of over-determining the CKM ma- 
trix to find signs of new physics [5] . The leptonic decay 
rate of the it via a W boson provides one such test. The 
QCD input to this is the pion decay constant, which is 
determined to 1% in lattice QCD [6]. If we take V u d from 
nuclear f3 decay [7] , we have a 2% determination of the 
leptonic decay rate to be compared to experiment. The 
leptonic decay rates of other charged pseudoscalars can 
also be determined to a few percent from lattice QCD [4J 
but then the comparison with experiment is generally 
needed to determine the appropriate CKM element. In- 
dependent tests of matrix elements, without CKM uncer- 
tainties, come only from electromagnetic decays. Here we 
provide two such tests through two different decay rates 
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of the J/ip: annihilation to e + e~ via a photon and ra- 
diative decay to the r\ c . We give the first results from 
full lattice QCD including u, d and s quarks in the sea, 
although earlier calculations have been done in quenched 
QCD [8] and including u and d sea quarks [10] . 

We are able to determine these matrix elements to a 
few percent because of our development of an accurate 
and fully relativistic approach to c quarks (as well as u, d 
and s) in lattice QCD called the Highly Improved Stag- 
gered Quark (HISQ) formalism [TT]. In this formalism we 
are able to normalise the vector current which mediates 
the electromagnetic decay accurately and nonperturba- 
tively and we show how to do that here. 

The layout of the paper is as follows: secti on [H] de- 
scribes the lattice calculation and then section |III| gives 
results for the J/ip mass, leptonic width (along with 
time moments of the J/ip correlator) and radiative de- 
cay rate in turn. We compare our results to experiment 
and to previous lattice QCD calculations in section |IV| 
Section |V| gives our conclusions. The Appendices discuss 
the more technical issues of discretisation errors and our 
two different methods for current renormalisation. 



II. LATTICE CALCULATION 

We use 6 ensembles of lattice gluon configurations at 4 
different, widely separated, values of the lattice spacing, 
provided by the MILC collaboration 12 . The configu- 
rations include the effect of u, d and s quarks in the sea 
with the improved staggered (asqtad) formalism. The u 
and d masses are taken to be the same with m u / d /m s 
approximately 0.2 on most of the ensembles. Based on 
our experience of other gold-plated mesons |2j we expect 
sea quark mass effects to be small for the J/ tp because it 
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Set ri/o 


auom l 


au m a a sq 


L s /a 


L t /a 


Sxi Sx s 


1 


2.647(3) 


0.005 


0.05 


24 


64 


0.11 0.43 


2 


2.618(3) 


0.01 


0.05 


20 


64 


0.25 0.43 


3 


2.658(3) 


0.01 


0.03 


20 


64 


0.25 -0.14 


4 


3.699(3) 


0.0062 


0.031 


28 


96 


0.20 0.19 


5 


5.296(7) 


0.0036 


0.018 


48 


144 


0.16 -0.03 


6 


7.115(20) 0.0028 


0.014 


64 


192 


0.17 0.04 



TABLE I: Ensembles (sets) of MILC configurations used for 
this analysis. The sea asqtad quark masses mf sq (I = u/d) 
and m"' are given in the MILC convention where Mo is the 
plaquette tadpole parameter. The lattice spacing values in 
units of ri after 'smoothing' are given in the second col- 
umn [12]. Here sets 1, 2 and 3 are 'coarse'; set 4, 'fine'; set 
5 'superfine' and set 6 'ultrafme'. The size of the lattices is 
given by L 3 S x L t . The final two columns give the difference 
between the sea quark mass and its physical value in units of 
the s quark mass [2]. 



has no valence light quarks. We can test this by compar- 
ison of results on sets 1 and 2 where the sea value of m u ^ 
changes by a factor of two and with set 3 where the sea 
value of m s changes by 70%. Table [i] lists the parameters 
of the ensembles. 

The lattice spacing is determined on an ensemble-by- 
ensemble basis using a parameter r\ that comes from fits 
to the static quark potential calculated on the lattice [H] . 
This parameter has small statistical/fitting errors but its 
physical value is not accessible to experiment and so must 
be determined using other quantities, calculated on the 
lattice, that are. We have determined r\ — 0.3133(23) fm 
using four different quantities ranging from the (2S-1S) 
splitting in the T system to the decay constant of the 
•q s (fixing fx and f n from experiment) [13] . Using our 
value for v\ and the MILC values for r\j a given in Tabled 
we can determine a in fm on each ensemble or, equiva- 
lently, a -1 in GeV needed to convert lattice quantities to 
physical units. 

On these ensembles we calculate c quark propagators 
using the HISQ action and combine them into meson 
correlation functions. The quark propagators are made 
from a 'random wall' source - a color 3-vector of U(l) 
random numbers - on a given timeslice to reduce the 
statistical noise. An added reduction comes from the 
use of a random starting point for the equally spaced 
time-sources we use on the coarse and fine ensembles. 
We include only connected correlation functions here - 
disconnected contributions for the J/ip are related to its 
hadronic width which is in keV and therefore negligible 
here. 

The c quark mass is tuned from the r\ c meson mass [2] . 
The appropriate 'experimental' mass for the r\ c for our 
calculations is 2.986(3) GeV, differing from the exper- 
imental result of 2.981(1) GeV [7] because of missing 
electromagnetic, r\ c annihilation and c-in-the-sea effects 
that we estimate perturbatively [14]. The HISQ lattice c 
quark masses for the ensembles we are using were deter- 
mined in [2]. 



Meson masses and decay constants are determined 
from simple '2-point' meson correlation functions made 
from combining quark propagators with appropriate spin 
matrices at source and sink to project onto the correct 
J PC . For staggered quarks, where the spin degree of free- 
dom has disappeared, the spin projection matrices are re- 
placed with space-time-dependent phases of ±1. Because 
of fermion-doubling, there are in fact 16 'tastes' of every 
meson made by combining a point-splitting of the quark 
and antiquark source and sink along with the appropriate 
±1 phases. The most accurate meson correlation func- 
tions come from either local or 1-link separated sources 
and sinks and we will restrict ourselves to these here. Be- 
cause the taste-splittings are discretisation effects we are 
free to use whichever taste is the most convenient for a 
given calculation. 

For the pseudoscalar mesons the mass differences be- 
tween the different tastes have a simple picture with the 
mass increasing as the amount of point-splitting in the 
source/sink operator increases. The lightest mass parti- 
cle is the Goldstone meson whose correlator is simply the 
modulus squared of the propagator and whose squared 
mass vanishes linearly with the quark mass. This is the 
one that is used to tune the quark mass. The other taste 
pseudoscalar mesons have a mass for which the difference 
of mass-squared with the Goldstone meson is a constant 
with quark mass which vanishes as a 2 s a 2 . These taste- 
splitting discretisation errors are particularly small with 
the HISQ action [TT]. They also become smaller, in pro- 
portion to the meson mass, as the meson mass increases 
and so are very small for mesons made of c quarks [llj . 
The mass difference between the Goldstone meson and 
the next heaviest pseudoscalar meson is visible, however. 
Both masses can be determined very accurately in lat- 
tice QCD because they both correspond to local opera- 
tors. The Goldstone meson corresponds to the local 75 
operator and the local non-Goldstone to the local 7075 
operator. We will use both of these mesons in our calcu- 
lation of the radiative decay rate of the J/ip. 

Vector meson taste-splittings are significantly smaller 
than for pseudoscalars and typically not visible for light 
mesons above the statistical errors. For the charmonium 
vectors we use the local 7, operator to determine the lep- 
tonic decay rate and two different 1-link split operators 
for the radiative decay. We discuss mass differences from 
taste-splittings further in Appendix [A"] 

III. RESULTS 

A. M J/4 , 

The determination of the mass of the J ftp is most ac- 
curately done through the determination of the charmo- 
nium hypcrfinc splitting, i.e. the mass difference with the 
pseudoscalar rj c meson. For the rj c we use the Goldstone 
meson, as discussed in section[Hj because this is the most 
accurately determined in lattice QCD and is the meson 
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we use to fix the c quark mass. We studied this meson in 
detail in [2]. For the J/i/j we use the local ji operator to 
create and destroy the vector meson. The J/ip correla- 
tors are then obtained by combining quark propagators 
from the default random wall with antiquark propagators 
from a source using the same random wall but patterned 
with phases, for example (— l) x for the vector polarised in 
the x direction. (— l) x is also inserted at the sink where 
the propagators are tied together. 

The J/ip and rj c correlators at zero spatial momentum 
are fit simultaneously so that correlations between them 
are taken into account. The fit form for the average J/ip 
correlator as a function of time separation between source 
and sink, t, is: 

C 2pt (t) = J2 a tHM tn ,t) - aijo{M io ,t) (1) 

with 

fn(M,i) = e -Mt + e -M(L t -t) 

fo(M,t) = (-l)*/ a fn(Af,t) (2) 

and L t the time extent of the lattice. i n — is the 
ground state and larger i n values denote radial or other 
excitations with the same J PC quantum numbers. The 
Mi n are the masses of the corresponding particles. There 
are 'oscillating' terms coming from opposite parity states, 
denoted i . The Goldstone i] c meson has the same fit 
form except that there are no oscillating contributions 
(when the r\ c is at rest). Note that we do not use any 
'smearing' functions for the propagators at either source 
or sink. 

To fit we use a number of exponentials i n , and where 
appropriate i , in the range 2-6, loosely constraining the 
higher order exponentials by the use of Bayesian pri- 
ors [15]. As the number of exponentials increases, we 
see the \ 2 value fall below 1 and the results for the fit- 
ted values and errors for the parameters for the ground 
state i = stabilise. This allows us to determine the 
ground state parameters ao and M as accurately as pos- 
sible whilst including the full systematic error from the 
presence of higher excitations in the correlation func- 
tion. We take the fit parameters to be the logarithm of 
the ground state masses M and M and the logarithms 
of the differences in mass between successive radial ex- 
citations (which are then forced to be positive). The 
Bayesian prior value for Mq for the r\ c is obtained from 
a simple 'effective mass' in the correlator and the prior 
width on the value is taken as 0.3. The prior value on 
M for the J/ip is taken to be 100 ± 50 MeV above the 
r) c . The prior value for mass splittings to and between 
excitations is taken as 600(300) MeV. The amplitudes 
a,i n and cii o are given prior widths of 1.0. We apply a cut 
on the range of eigenvalues from the correlation matrix 
that are used in the fit of 10~ 4 . We also cut out small 
t/a (and (L t —t)/a) values below 6 from our fit to reduce 
the effect of higher excitations. 



*V 0.027 

\ 0.026 
+ 

J 0.025 
% 0.024 
j 0.023 





1 1 1 


1 1 1 1 1 




1 


> 














1 1 1 1 1 1 1 1 1 



10 20 30 40 50 60 70 80 90 

t/a 



FIG. 1: Our average J/f/j correlator divided by the ground 
state exponential (fn(Mo, i) from eq. |2|) as a function of lat- 
tice time. Lines are drawn to join the points (which include 
statistical errors) for clarity. The fitted result for the ground 
state amplitude, a$, is given by the blue band. The fit in- 
cludes 6 normal exponentials and 6 oscillating ones, which 
are responsible for the oscillating behaviour clearly seen in 
the results. 



Figure [T] shows the quality of our results with a plot of 
the J/ip correlation function. It is divided by the ground- 
state exponential function so that it shows a plateau in 
the centre of value Oq. The results for the ground-state 
masses in lattice units of the J jt\> and r\ c and the dif- 
ference between them, aAMh yp , are given in Table [TT1 
The difference is typically more accurate than that ot> 
tained by simply subtracting the masses because of the 
correlation between the correlators. 

The hyperfine splitting is converted to physical units 
using the values for a on each ensemble as discussed in 
section [TTJ. The results are shown in Figure [2j Figure [2] 
includes the error from the determination of the lattice 
spacing on each point. This dominates the error but is 
correlated between the points and that should be borne 
in mind in looking at the figure. It is important to re- 
alise that the naive lattice spacing error is magnified by a 
factor of approximately two in the hyperfine splitting be- 
cause of the inverse relationship between hyperfine split- 
ting and quark mass. For example, a shift by uncertainty 
S upwards in the inverse lattice spacing causes a shift 
upwards in the meson mass by the same proportion. To 
determine the total effect of this on the hyperfine split- 
ting we must include the effect of retuning the c quark 
mass to make the meson mass correct again. This means 
in this case retuning the quark mass down by fraction 5 
which shifts the hyperfine splitting upward by a further 
factor of S to that coming simply from the lattice spacing 
change. Thus the change in the hyperfine splitting, rep- 
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Set TVcfg x N t 


m c a 


e aM Vc 


aMj/nf, aAMhyp afj/^/Z Z cc 


1 


2099 x 8 


0.622 


-0.221 1.79118(4) 


1.85934(8) 0.06817(6) 0.2810(2) 0.979(12) 


2 
2 


2259 x 4 
2259 x 8 


0.63 

U.DD 


-0.226 1.80851(5) 
-0.244 1.86667(4) 


1.87797(10) 0.06946(8) 0.2855(2) 0.979(12) 
1.93430(9) 0.06763(7) 0.2925(2) 0.974(12) 


3 


323 x 8 


0.617 


-0.218 1.78212(12) 1.85081(23) 0.06869(17) 0.2804(5) 0.979(12) 


4 


566 x 4 


0.413 


-0.107 1.28052(7) 


1.32901(12) 0.04849(10) 0.1829(2) 0.983(12) 


5 


200 x 2 


0.273 


-0.0487 0.89948(8) 


0.93369(13) 0.03421(11) 0.1244(3) 0.986(12) 


6 


208 x 1 


0.193 


-0.0247 0.66649(6) 


0.69217(11) 0.02568(10) 0.0925(3) 0.990(12) 



TABLE II: Results in lattice units for the masses of r\ c and J/ip and their difference on each ensemble along with the raw 
(unrenormalised) decay constant and Z factor for the J/ip. Columns 3 and 4 give the bare HISQ charm quark mass, tuned 
from the rj c and the corresponding coefficient e used in the Naik discretization improvement term of the HISQ action 2 . All of 
the charm quark masses are very well tuned except for the lower result on set 2 (m c a — 0.66), which was deliberately mistuned 
to assess the sensitivity of quantities to the tuning. Of the remaining masses the least well-tuned is on superfine set 5 where 
M 7)c is 0.5% too high. Column 2 gives the number of configurations used and the number of time sources for propagators on 
each configuration. Results are binned on time sources and binned over neighbouring configurations for sets 5 and 6. The J/ip 
correlators are averaged over polarisations except on sets 2 and 3 where only one polarisation was calculated. The results for 
the r/ c masses are also given in [2]. They differ slightly from these in some cases because of fitting simultaneously with J/tp 
correlators. The Z factors are taken from moment 4 of the nonperturbative (on the lattice) current-current correlator method 
described in Appendix |B 1| 
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FIG. 2: Results for the charmonium hyperfine splitting plot- 
ted as a function of lattice spacing. For the a;-axis we use 
(m c a) 2 to allow the a-dependence of our fit function (eq. (JsJ) ) 
(blue dashed line with grey error band) to be displayed sim- 
ply. The data points have been corrected for c quark mass 
mistuning and sea quark mass effects, but the corrections are 
smaller than the error bars. We do not include on the plot 
the deliberately mistuned c mass but it is included in the 
fit to constrain the c mass dependence. The errors shown 
include (and are dominated by) uncertainties from the de- 
termination of the lattice spacing a (from the physical value 
of the parameter ri) that are correlated between the points. 
The experimental average is plotted as the black point at the 
origin, offset slightly from the y-axis for clarity. 



resenting its uncertainty, is approximately 28 |16j 1 . Thus 
lattice spacing uncertainties are typically much more im- 
portant in the determination of hyperfine splittings than 



1 This point has frequently been overlooked in lattice QCD calcu- 
lations. 



statistical errors. 

We fit the hyperfine splitting as a function of lattice 
spacing and sea quark masses to the form: 

f(a,8xi,5x s ) = f x (3) 
2^Ciiw(on*c) (^y(^) (^) 

ijkl 

+ (d Q + di(am c ) 2 )(M VcM tt - M Vctexpt ). 



Here fo is the physical result, the sum over ijkl allows for 
discretisation errors and sea quark effects and the final 
term allows for mistuning of the c quark mass. We allow 
the discretisation errors, which are evident in our results, 
to have a scale set by the c quark mass. These appear 
only as even powers of a for staggered quarks. Sxi and 
8x s are the mistuning of the sea quark masses: 



^s,phys 



q,phys 



(4) 



Sxi and Sx s values are given for each ensemble in Ta- 
ble [T] and are taken from Appendix A of [2]. Eq. ^ 
includes a term for each sea quark [u/d appearing twice, 
and s), with the coefficients constrained to be the same 
so that the fit function is symmetric with respect to inter- 
change of any two. The division by 10 is because the scale 
for dependence on light quark masses from chiral pertur- 
bation theory is 4^/^ sw \Qm s . We sec no significant 
sea quark mass dependence in the hyperfine splitting. A 
fairly strong dependence was seen in the twisted mass cal- 
culations [TU]. However, at least some of that dependence 
could be attributed to the sea quark mass dependence of 
the lattice spacing, since that is determined only in the 
chiral limit. Here we determine the lattice spacing for 
each ensemble and hence separate lattice spacing depen- 
dence from physical sea quark mass effects. The sum over 
ijkl in eq. ([3| allows for the possibility of lattice spacing 
dependent sea quark mass effects. 
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We take a Bayesian prior [T5] on /o of 0.1(1) and then 
fix Coooo to 1. The other cyfcj are given priors of 0.0 ± 1.0 
except for the Cojki which determine the a-independcnt 
sea quark mass dependence. These are taken to have pri- 
ors 0.0±0.33 because we expect sea quark mass effects to 
be typically a factor of 3 smaller than valence quark mass 
effects which would have chiral perturbation theory coef- 
ficients of 0(1). We include 5 terms in the a-dependence 
and 3 in the 5x dependence. Including additional terms 
makes no difference to the value for /o or its error. The 
priors for d an d di are taken as 0.00(5), informed by 
the expectation that the hyperfine splitting should be in- 
versely proportional to the mass, and by the effect of our 
mistuned c mass on set 2 which agrees roughly with that 
expectation. 

The fit gives f = 116.5(2. l)MeV, as the result for the 
hyperfine splitting in the absence of clcctromagnetism, c- 
in-thc-sca and cc annihilation. The first two affect the rj c 
and J/tjj equally and so have no effect on the hyperfine 
splitting. The third affects the rj c more than the J/ip, 
which has negligible width. A perturbative estimate of 
the shift of the rj c mass resulting from its annihilation 
to two gluons [TT] related this to the total i] c width and 
obtained a shift downwards of the r\ c mass of 2.4 McV 2 . 
Using this, we have since applied a shift of 2.4(1.2) MeV 
for this effect to determine the rj c mass to which to tune 
our c quark mass, as in section [IT] For this purpose the 
impact of the shift is completely negligible, amounting to 
less than 0.1% of the r\ c mass. For the hyperfine split- 
ting, however, this shift could be a relatively large effect. 
Nonperturbative calculations of the contribution of 'dis- 
connected diagrams' to the r\ c mass have agreed on a 
small value of a few MeV for the shift from r] c annihila- 
tion but obtained the opposite sign [17]. The argument 
is that the perturbative result may be modified signifi- 
cantly by the gg intermediate state forming a resonance 
such as a glueball which is lighter in mass than the r) c , or 
a lighter hadron state. To allow for this possibility and 
be consistent with the nonperturbative calculations we 
do not apply a shift to the hyperfine splitting obtained 
from our fit above, but instead take an additional sys- 
tematic error of 2.4 MeV, corresponding to our original 
shift, to allow for the effect. 

Our final result for the hyperfine splitting is then: 

AM hyp = 116.5(2.1) (2.4) MeV (5) 

where the errors are in turn from statistics/fitting and 
t] c annihilation. The uncertainty from rj c annihilation 
dominates the error. A complete error budget is given in 
Table HTU 

This is to be compared to the difference of the ex- 
perimental averages of the two masses of 115.9(1.1) 



2 This would now amount to 2.9 MeV given that the average ex- 
perimental width of the rj c has increased to 30 MeV [7]. 



Mj/if, — M Vt , fj/if, Vj/if,^ Val (0) 



(am c ) z extrapolation 


0.45 


0.45 


3.5 


statistics 


0.50 


0.41 


0.74 


lattice spacing 


1.6 


0.42 


0.0 


sea quark extrapolation 


0.29 


0.26 


1.3 


M Vl . tuning 


0.11 


0.09 


0.0 


Z 




1.23 


0.14 


M Vc annihilation 


2.1 


0.0 


0.0 


electromagnetism 


0.0 


0.5 


0.5 


Total (%) 


2.7 


1.5 


3.8 



TABLE III: Complete error budget for hyperfine splitting, 
leptonic width and vector form factor as a percentage of the 
final answer. 

MeV [7J. Quite a spread of results make up the av- 
erage. Recent values tend to be at the lower end of 
the hyperfine splitting range. For example, the 2011 
Belle result for the rj c mass gives a hyperfine splitting 
of 111.5(+^) MeV Q2], and a recent result from BESIII 
gives 112.6(0.9) MeV [15]. 

B. r(J/V> -» e+e") and i? e+e - 

The amplitude, ao, from the fit in equation ([I]) to our 
J/ip correlators is directly related to the matrix element 
for the local vector operator to create or destroy the 
ground-state vector meson from the vacuum. The vector 
meson decay constant, /„, for meson v is defined by: 

(0|^yVl«> = fvmvf (6) 

where e l is the meson polarization. f v for the J/tp is then 
determined from our lattice QCD correlators, in terms of 
the ground-state parameters from our fit (eq. ([T])) by: 




where Z is the renormalisation constant required to 
match the local vector current in lattice QCD used here 
to that of continuum QCD at each value of the lattice 
spacing. 

f v is clearly a measure of the internal structure of a 
meson and in turn is related to the experimentally mea- 
surable leptonic branching fraction: 

r(„^eV)4«W- ( g ) 

where is the electric charge of the heavy quark in units 
of e (2/3 for c). The experimental average, T(J/ip — > 
e+e") = 5.55(14)keV [7] gives / J/v , = 407(5) MeV, re- 
membering that the electromagnetic coupling constant 
runs with scale and using \/aQEDijn c ) = 134 [20] , This 
can then provide a test of QCD at the 1% level. Electro- 
magnetic corrections are small since the J/ip must decay 
to an odd number of photons [21] . 
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Our results for fj/^/Z are given in Table \U\ The fi- 
nal column of that table gives the values of Z deter- 
mined from current-current correlators as described in 
Appendix [B] This method uses continuum perturbation 
theory through 0(a^) to normalise the lattice QCD cor- 
relators at small times. Z then results from a combina- 
tion of non-perturbative lattice QCD calculations with 
continuum perturbation theory in a similar approach to 
that of the RI-MOM scheme 3 used to renormalise the cur- 
rents for the same calculation using twisted mass quarks 
in [10] . The current-current correlator method has the 
advantage that we can use the same correlators from 
which we also extract, at large times, the nonpertur- 
bative information on the ground-state mass and decay 
constant. Indeed this allows some cancellation of dis- 
cretisation errors apparent in the unrenormalized decay 
constant. 

Multiplying fj/.^/Z by Z and then by a" 1 in GeV 
gives the physical results for the decay constant plotted 
in Figure [3] We fit these to the same function of lat- 
tice spacing and sea quark mass used for the hypcrfine 
splitting, eq. ([3]). The only differences are that the prior 
on /o is taken as 0.5(5) in this case and the priors on 
the slope of the variation of fj/^p with M 7]c are taken as: 
d , 0.065(5) and d u 0.00(25). These are informed by the 
variation we see for the deliberately mistuned c mass on 
set 2 and also by our extensive study of the behaviour of 
f ric with M 7)c in [5] . There we find a strong a-dependence 
in the slope of the decay constant with mass and so we 
allow for that here. 

The physical result that we obtain in the continuum 
limit is: 



f m = 405(6) (2)MeV. 



(9) 



The first error is from the fit and is dominated by the 
error from the Z factor. The second error is an estimate 
of systematic effects from missing electromagnetism in 
our lattice QCD calculation [2 . The effect of missing 
c-in-the-sea is negligible in this case. A complete error 
budget is given in Table |III| 

The leptonic width is determined by the amplitude of 
the ground-state that dominates the correlator at large 
times. We can also determine the charm contribution 
to R e + e - through the time moments of the J/ip correla- 
tor which depend on the behaviour at short times. The 
moments are defined by: 



2 E" 



(10) 



where t is lattice time symmetrised around the centre of 
the lattice (see Appendix^. Results for (G^/Z 2 ) 1/(n ~ 2) 
in lattice units on each of our ensembles are given in Ta- 
ble IW] for n = 4,6,8 and 10. The power l/(n - 2) is 




FIG. 3: Results for the charmonium vector decay constant 
plotted as a function of lattice spacing. For the x-axis we use 
(m c a) 2 to allow the a-dependence of our fit function (eq. |3j)) 
(blue dashed line with grey error band) to be displayed sim- 
ply. The data points have been corrected for c quark mass 
mistiming and sea quark mass effects, but the corrections are 
smaller than the error bars. We do not include on the plot the 
deliberately mistuned c mass but it is included in the fit to 
constrain the c mass dependence. The errors shown include 
(and are dominated by) uncertainties from the determination 
of the current renormalization factor, Z, that are correlated 
between the points. The experimental average is plotted as 
the black point at the origin, offset slightly from the y-axis 
for clarity. 



taken to reduce all the moments to the same dimension. 
We take the Z factor for the vector current to be the 
same one used for the leptonic width above, determined 
in Appendix [B] Figure [4] then shows the physical results 
for the moments as a function of lattice spacing. The 
gray bands show our fits which use the same function of 
lattice spacing and sea quark masses as given in eq. (J3j) . 
We reduce the prior width on the lattice spacing depen- 
dent terms by a factor of 4 because the moments are not 
as sensitive to short distances as the leptonic width or 
hyperfine splitting. 

The physical results that we obtain for each moment 
in the continuum limit are given by: 



(gD 1/2 

(G 6 y ) 1/4 

(gD 1/6 

(cr ) i/8 



0.3152(41)(9) GeV" 1 
0.6695(57)(13) GeV" 1 
0.9967(65)(10) GeV" 1 
1.3050(65)(6)GeV _1 . 



(11) 



3 This method is often called 'nonperturbative' in the lattice QCD 
literature. 



The first error comes from the fit and the second allows 
for electromagnetism (e.g. photons in the final state) 
missing from our calculation but present in experiment. 
The error is estimated by substituting olqed fo r ol s in 
the perturbative QCD analysis of the moments [55]. A 
complete error budget for our results is given in Table |V| 
The results agree well with the values extracted for the 
q 2 derivative moments, M.]~ (n = 2k + 2), of the charm 
quark vacuum polarization using experimental values for 
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Set 


/ gY \ 1/2 
m c a \-^) 


( gy y/ 4 


1 gY \ 1/6 

\ Z^a® ) 


( gYo \ 1/8 

\Z^a» ) 


1 


0.622 0.5399(1) 


1.2162(1) 


1.7732(1) 


2.2780(1) 


2 


0.63 0.5339(1) 


1.2054(1) 


1.7581(1) 


2.2584(1) 


2 


0.66 0.5135(1) 


l.l692(l) 


l.708l(l) 


2.1941(1) 


3 


0.617 0.5434(1) 


l. 2223(1) 


1.7817(1) 


2.2888(1) 


4 


0.413 0.7586(1) 


1.6351(1) 


2.3887(2) 


3.0952(2) 


5 


0.273 1.0681(1) 


2.2705(2) 


3.3454(3) 


4.3601(4) 


6 


0.193 1.4323(3) 


3.0397(5) 


4.4990(7) 


5.8738(8) 



TABLE IV: Results in lattice units for time moments of the 
J ftp correlator as defined in eq. < |10[ ). We give results for n=4, 
6, 8 and 10. 





(Gl) 1/2 


(G 6 T /4 


(G 8 Y /6 


(cr ) i/8 


(am c ) 2 extrapolation 


0.18 


0.18 


0.16 


0.16 


statistics 


0.05 


0.04 


0.03 


0.03 


lattice spacing 


0.32 


0.51 


0.43 


0.30 


sea quark extrapolation 


0.14 


0.13 


0.12 


0.12 


M Vo tuning 


0.15 


0.18 


0.17 


0.16 


Z 


1.23 


0.61 


0.41 


0.31 


electromagnetism 


0.3 


0.2 


0.1 


0.05 


Total (%) 


1.3 


0.9 


0.7 


0.5 



TABLE V: Complete error budget for the time moments of 
the J/ip correlator as a percentage of the final answer. 



R e + e - = <r(e + e~ — > haclrons)/cr pt (221 [23]. The values, 
extracted from experiment by |22| and appropriately nor- 
malised for the comparison to ours, are: 

(M 1 cxp 4!/(127r 2 e 2)) 1 / 2 = 0.3142(22) GeV" 1 

(M 2 cxp 6!/(127r 2 e 2 )) 1/4 = 0.6727(30) GeV" 1 

(M 3 oxp 8!/(127r 2 e 2 )) 1/6 = 1.0008(34) GeV" 1 

(Mf p 10!/(127r 2 e 2 )) 1/8 = 1.3088(35) GeV _1 . (12) 

Our results from lattice QCD have approximately double 
the error of the experimental values but together these 
results provide a further test of QCD to better than 1.5%. 



c. r(j/v> -+ jvc) 

The radiative decay of the J/tp meson to the r\ c re- 
quires the emission of a photon from either the charm 
quark or antiquark and a spin-flip, so it is an Ml transi- 
tion. Because it is sensitive to relativistic corrections this 
rate is hard to predict in nonrelativistic effective theories 
and potential models (see, for example, [2~4"l |2"5] ) Here 
we use a fully relativistic method in lattice QCD with 
a nonpcrturbatively determined current renormalisation 
and so none of these issues apply. In addition, of course, 
the lattice QCD result is free from model-dependence. 

The quantity that parameterises the nonperturbative 
QCD information (akin to the decay constant of the pre- 
vious section) is the vector form factor, V(q 2 ), where q 2 
is the square of the 4-momentum transfer from J/if> to 



1.4 



^ 1-2 

i 

> 

A 0.8 



I 0.6 
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£ 

Is 0- 4 
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n = 6 




• 
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1 


n = 4 
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FIG. 4: Results for the 4th, 6th, 8th and 10th time moments 
of the charmonium vector correlator shown as blue points and 
plotted as a function of lattice spacing. The errors shown (the 
same size or smaller than the points) include (and are domi- 
nated by) uncertainties from the determination of the current 
renormalization factor, Z, that are correlated between the 
points. The data points have been corrected for c quark mass 
mistuning and sea quark mass effects, but the corrections are 
smaller than the error bars (the value for the deliberately 
mistuned c mass on set 2 is not shown). The blue dashed 
line with grey error band displays our continuum/chiral fit. 
Experimental results determined from R e + e - (eq. ( |12[ )) are 
plotted as the black points at the origin offset slightly from 
the y-axis for clarity. 



rj c . The form factor is related to the matrix element of 
the vector current between the two mesons by: 



{Vc(p')\crfc\J/^(p)) = 



2V(q 2 ) 



(13) 

Note that the right-hand-side vanishes unless all the vec- 
tors are in different directions. Here we use a normalisa- 
tion for V(q 2 ) appropriate to a lattice QCD calculation 
in which the vector current is inserted in one c quark line 
only and the quark electric charge (2e/3) is taken as a 
separate factor. The decay rate is then given by [8]: 



r(J/V> -> 77 c 7) = Q.QED 



27(M„ c + Mj^y 



: \V(0)\ 2 , (14) 



where it is the form factor at q 2 = that contributes be- 
cause the real photon is massless. |<j| is the corresponding 
momentum of the rj c in the J/tp rest-frame. 
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t 



Vc 




V 




J/i> 



T 



FIG. 5: A schematic diagram of the connected '3-point' func- 
tion in lattice QCD for J/ip to r) c radiative decay. The lines 
all represent c quark propagators in this case. The propagator 
labelled 1 is the spectator quark; 2 and 3 are the initial and 
final active quarks respectively. and T label the position 
in time of the r\ c and J/ip operators. The vector current is 
inserted at time t which takes all values between and T. 



The most recent experimental result from CLEO-c [26] 
of 1.98(31)% for the branching fraction, combined with 
the total width of the J/ip of 92.9(2.8) keV [7] gives 



V(0) 



GXpt 



1.63(14), 



(15) 



where we have used ozqed — 1/137 and |g| = {Mj/^ — 
M Vc )(M J / i> + M Vc )/(2Mj/^), The value of \q\ from ex- 
periment is 0.1137(11) GeV where the error comes from 
the uncertainty in the r\ c mass. V(0) is then the quantity 
that can be calculated in lattice QCD and compared to 
experiment. 

The radiative decay of the J/ijj to rj c meson needs the 
calculation of a '3-point' function in lattice QCD. The 
3 points (in lattice time) correspond to: the position of 
the rj c operator, which we take as the origin; the position 
of the J I ip operator which we denote T and the position 
of the insertion of a vector operator, V = C7 M c, which 
couples to the photon at time t. t varies from to T. 
Sums over spatial points are implied at each time. The 
'connected' correlator that we calculate is illustrated in 
Figure [5] Disconnected correlators are expected to be 
negligible here based on perturbative and phenomeno- 
logical arguments [8] and we do not include them. 

The 3-point function is calculated in lattice QCD by 
combining 3 quark propagators together with appropri- 
ate spin projection matrices. As discussed in section [n] 
for staggered quarks these 7 matrices become ±1 phases. 
Tastes must be combined in a staggered quark correla- 
tor so that the overall correlation function is 'tasteless'. 
What this means for a 3-point function is that only cer- 
tain taste combinations of J/ip, Vc and V operators are 
allowed. To optimise statistical errors we need to keep to 
a minimum the amount of point-splitting in the opera- 
tors. It is also convenient, for renormalisation purposes, 
to have a vector current, V, which corresponds to a lo- 
cal operator (and this is also what we used for the decay 
constant in section MB). 



We therefore choose the rj c operator to be the local 75 
operator (so that the r] c is the Goldstone pseudoscalar 
with spin-taste 75 <g> 75) and the J/ip operator to be a 
one-link separated 7o7i operator in which the polarisa- 
tion of the J/ip and the one-link separation are both in 
an orthogonal spatial direction to the polarisation of the 
vector current, V = crfkC (this J/ip has spin-taste struc- 
ture 70 7* ® 7o7i7j)- 

To implement this configuration is simple. The specta- 
tor quark propagator (number 1 in Figure pi is generated 
from the default random wall at time 0. Active propa- 
gator 2 is then generated from a source which is made 
from a symmetric point-splitting of propagator 1 at time 
T patterned by a phase. For a J/ip with polarisation 
x we take a point-splitting in the y direction and phase 
( — l) x+z . Active propagator 3 is made from the same 
default random wall as 1. Finally 2 and 3 are combined 
together at t by summing over space with a patterning 
of (— l) z to implement a local vector current in the z 
direction. 

To achieve the configuration corresponding to q 2 = 
we keep the J/ip at rest in the frame of the lattice and 
give the r\ c an appropriate spatial momentum. The rj c 
momentum is implemented by calculating propagator 3 
with a 'twisted boundary condition' |27U28| . If propaga- 
tor 3 is calculated with boundary condition: 



„2iri6. 



(16) 



then the momentum of the r\ c meson made by combining 
propagators 1 and 3 with our random wall sources and 
summing over spatial sites at the sink is: 



Pj 



2n o 



(17) 



The boundary condition in eq. (|16j) is actually imple- 
mented by multiplying the gluon links in the j direction 
by phase exp(2wi9j/L s ). We take j to be the y direction 
here so that the momentum is in an orthogonal direction 
to the polarisation of both the J/ip and V. 
The 3-point function is then given by: 

C 3pt (0,t,T) = \(-l) XT+ZT (-l) Zt x (18) 

S T ,S t 

Tr [g(t, T) [g(T + l y , 0) + g(T - l y , 0)]ffj(t, 0) } 

where g represent staggered c quark propagators, with 
gg computed with a phase on the gluon field, the trace 
is over color and sums are done over spatial sites St and 
St at t and T. The 1/4 is the taste factor for the nor- 
malisation of a staggered quark loop. The corresponding 
2-point function for the rj c meson is 



C Ve ,2pt(p,t) = 2 x Tr k*>°)4M) 



(19) 
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Set Ndg x N t 


T values 


m c a 


aM J/4 , 




2tt(9 aE u Va V n n V(0)/Z Z ff 


'2 2 

a q 










7o 7i ® 7o7i7j 


75 ® 75 






1 


2088 x 4 


15,18,21 


0.622 




1 7Q1 1 (Sl'4'l 


1 R4in 1 7Q94' : i('4 x i ()' : {p,'?('?\ 1 QOflHl"; flQRQfSl'11 > ) 


1 (A\ x 1 n~ b 


2 


2259 x 4 


15,18,21 
15,18,21 


0.63 
0.63 


1.87962(14) 


1.80839(8) 


1 4007 1 81fl9^(' c i x i n nw»f9l 1 8Q7l'19 s l fl Q8Q4('S , 1 

1.3880 1.81019(4) 0.0362(4) 1.883(20) 0.9894(8) 


_7Ci v 1 n~ b 
1(5) x 10~ 5 


4 


ifli i w A 

iyii x 4 


20,23,26,29 0.413 


1 ^9Q0 c il'Q > l 






U^4: J A 1U 










;i <> /» n 


[0 [5 W /U Jo 






1 


2088 x 4 


15,18,21 


0.622 


1.86035(15) 


1.79621(4) 


1.5120 1.79725(4) 0.0338(6) 1.925(35) 0.9896(11) 


3(5) x 10" 5 


2 


2259 x 4 


15,18,21 
15,18,21 


0.63 
0.66 


1.87887(13) 
1.93604(15) 


1.81369(6) 
1.87254(6) 


1.2814 1.81480(5) 0.0334(8) 1.896(45) 0.9894(8) 
1.2490 1.87355(6) 0.0322(8) 1.934(42) 0.9863(17) 


0(4) x 10~ b 
1(5) x 10" 5 


4 


1911 x 4 


19,20,23,26 0.413 


1.32904(11) 


1.28160(4) 


1.3116 1.28243(4) 0.0342(4) 1.872(21) 1.0049(10) 


-2(1) x 10" b 



TABLE VI: Results from simultaneous fits for 3-point and 2-point correlators for J/lp — > yr) c decay. The upper table gives 
results from our preferred jpsigammaO method; the lower table from etacgammaO. See the text for a definition of the two 
methods. Column 2 gives the number of configurations and time sources for on each configuration. Column 3 gives the 
different values of the end-point of the 3-point function, T, included in the fit. The lattice c quark mass and e parameter 



are the same as those used in section III A (the lower table includes the deliberately mistuned mass on set 2 for comparison). 
clMj/t), and aM Vc are the zero-momentum meson masses for the tastes of J/ip and rj c mesons used here. 2n8 indicates the 
value of the phase at the boundary used to achieve the kinematics of q 2 = in the J/ip — > r\ c decay. The a 2 q 2 values actually 
obtained with those kinematics are given in the final column (rows 2 and 3 of the upper table compare two different values of 
a 2 q 2 close to zer o). aE rfa gives the energy of the r\ c at the value of the spatial momentum corresponding to 0. V°° from the 
3-point fit of eq (251 is given in column 9 and this is converted to a value of V(0)/Z in column 10 using eq. ( 27[) . C olumn 11 
gives the values of the renormalisation parameter, Z, obtained from the vector form factor method of Appendix |B 2| 



The 2-point function for the J/ip is given by 



(20) 



Tr [g(t, O)0?t(i + l y , l y ) + g\t - l y , l y ) + {1 o -1})] 

As an alternative configuration we can take the r\ c op- 
erator to be the local 7075 operator (so that the r\ c is 
the local non-Goldstone meson with spin-taste structure 
7o75 ®7o75) and the J/ip operator to be a one-link sepa- 
rated 7j operator in which the polarisation of the J/ip and 
the one-link separation are both in an orthogonal spatial 
direction to the polarisation of the vector current, V (this 
has spin-taste structure ji ®7i7j). The 3-point function 
is then given by: 

C 3pt (0,i,T) = ]T I(-l)»o+«»+*(_l)«r(_l)*. x 

ST;S t 

Tr [g(t, T)(g(T + l y , 0) + g(T - l tf> 0))ffj(i, 0)] (21) 
and the corresponding 2-point functions are: 

C Vc ,2 P t{0,t) = J2 1 (- l) Xo+Vo+Zo (- l) Xt +vt +Zt x 



St 

Tr 



[g(t,0)gl(t,0) 



(22) 



and 



C 



J/4>,2pt 



(0,t) 



~y ' \^_-^x a +z +t ^_-^x t +z t +t t x ^23) 



Tr [g(t, 0)tf(t + l y , l y ) + g\ t - ly, ly) + {1 O -1})] 



We call this configuration the 'etacgammaO' configura- 
tion and the original configuration of eq. ( 19 1 the 'jp- 
sigammaO' configuration. In fact, as we shall see, the 
jpsigammaO configuration is to be preferred on the basis 
of statistical errors but the results agree between the two. 



The 3-point function in both cases is calculated along 
with the 2-point functions for the r\ c and J/ip mesons 
that appear in it. We use multiple time-sources for point 
on each configuration and also multiple values for T. 
Figure [6] shows results for the 3-point function on fine 
set 4, normalising it to the product of the relevant 2- 
point functions. The two plots compare results for the 
jpsigammaO method and the etacgammaO method. The 
two differ in the amount of oscillation that is seen at the 
two ends of the plot. Not surprisingly the jpsigammaO 
method shows more oscillation on the J/ip end (t near 
T) since the rj c in this case would not oscillate at rest. 
The etacgammaO method has relatively large oscillations 
for the r\ c side but smaller oscillations on the J ftp side. 
In both cases statistical errors are very small enabling 
us to fit both normal and oscillating terms. The figure 
also shows how having multiple T values improves our 
determination of the ground-state transition amplitude. 



We fit the 3-point function and 2-point functions si- 
multaneously to a multi-exponential that determines the 
ground-state amplitudes accurately because it includes 
excited state contributions. The fit form for the 2-point 
functions was already given in eq. |l]). Here both the J/ip 
and rj c correlators have oscillating contributions and, in 
the T] c case, the exponent gives the energy of the meson 
at momentum pj (eq. (17)) rather than the mass. The 
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and, again: 



FIG. 6: AplotoftheratioC , 3pt(t,T)/(C2 P t,, c (t)C 2p t,j/v.(T- 
t)) for the three values of T on the fine ensemble, set 4 and 
for our two different methods. Lines join the points (which 
have statistical errors on them) for clarity. We only include 
points in the central region of t i.e. t>5 or t<T — 5. The 
pink shaded band shows the ratio of fit parameters VooVeso&o 
which is the ground-state contribution to this ratio. These 
come from a fit that included 6 normal exponentials and 5 
oscillating ones (which produce the oscillations evident in the 
figure). The top plot shows the results for the case where 70 is 
included in the J/tj) operator (jpsigammaO method) and the 
lower plot shows the results for the case where 70 is included 
in the r/ c operator (etacgammaO method). 



fit form for the 3-point function is then: 

C 3pt (t,T)= (24) 
J2 air, HE a ,in > t)Vg\j n b jn in(E bdn ,T-t) 

- ^ME a ^,t)V™ Jo b 3 MEb,3 ,T-t) 

in, jo 

- J2 <HKio,t)V°X b 3MEb,j n ,T-t) 

io,3n 

+ <HEa,io,t)V^ jo b j MEb, jo ,T - t) 



in(E,t) = e- Et + e- E ^-^ 
io(E,t) = (-l) t/a in(E,t) 



(25) 



with L t again the time extent of the lattice. Here n 
denotes the normal contributions and o the contribu- 
tions from oscillating states. The ground-state ener- 
gies/masses that we need are £"^ C) o and Eji^ t0 = Mj/^ 
and the matrix element between them that is propor- 
tional to Vo"o • By matching to a continuum correlator 
with a relativistic normalisation of states and allowing 
for a renormalisation of the lattice vector current we see 
that: 



( Vc \V\Jm = 2Z^ 



M J/-4> E Vc V Qfi- 



(26) 



The vector form factor that we need, V(0), is then, from 
eq. (Il3| , given by: 



no) 

z 



M 



2M 



Mr, 
J/ipPj 



Mj^E v V™, (27) 



Wo 



with pj from eq. (17 1. The determination of Z will be 
discussed below. 

To perform the joint fit to the 3pt correlators using 
eq. (25) and the 2pt correlators using eq.(T]) we use the 
same approach as outlined in section |III A For both the 
r\ c and J/ip, the prior for the ground-state mass comes 
from the effective mass of the correlator. We use priors 
of 600 MeV with a width of 300 MeV for the difference 
in mass between the ground state and the lowest oscillat- 
ing mass and between all radial excitations, both normal 
and oscillating. The 2-point amplitudes, ai and bi, have 
prior widths of 0.5 and the 3-point amplitudes, V^, have 
widths of 0.25. We omit t values below a certain t m ; n to 
reduce the effect of excited states. t mln = 4(5) for the 
coarse(fine) lattices for the etacgammO method and 6 for 
the jpsigammO method. 

Table |VI| gives our results from fits that include 6 nor- 
mal exponentials and 5 oscillating. We work on ensem- 
bles 1, 2 and 4 of Table [I] but using more configurations 
than in section llH Al to reduce statistical errors. Table fVTl 
gives the number of configurations and time sources as 
well as the values of T used in the 3-point functions. 
It is important to use both even and odd values of T to 
separate clearly the normal and oscillating contributions. 
Having determined the mass of the local non-Goldstone 
rj c and 1-link vector from separate 2-point function fits we 
then determine the value of 9 needed to achieve q 2 = 0. 
The final fits are done as a simultaneous fit to the 3-point 
function and 2-point functions for zero momentum and 
finite momentum rj c and zero momentum J/ip. 

The key parameters to be determined from the fit, as 
discussed above, are the ground-state masses of the r\ c 
and J/ip, the ground state energy at non-zero momen- 
tum of the 7/ c and the ground-state to ground-state am- 
plitude of the 3-point function. Our fit returns excited 
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FIG. 7: Results for the vector form factor at q 2 = for J/^ — > 
r/ c decay plotted as a function of lattice spacing. The filled 
blue circles are from our preferred jpsigammaO method; the 
open blue circles are from the etacgammaO method. For the 
a;-axis we use (m c a) 2 to allow the a-dependence of our fit 
function to be displayed simply (blue dashed line and grey 
band). The fit is to results from the jpsigammaO method. The 
errors shown include statistical errors and errors from the Z 
factor. The experimental result extracted from the branching 
fraction for J/tp — > ^r\ c is plotted as the black point offset 
slightly from the origin for clarity. 



state to ground-state and oscillating to ground state am- 
plitudes also. Most of these do not have a significant 
signal. Indeed the excited state to ground state ampli- 
tudes are very small, as expected since they correspond 
to a hindered Ml transition. A non-zero result is seen 
for the transition between the oscillating partner of the 
r] c (in the etacgammaO method) and the J/ip- This corre- 
sponds to the El Xco to J/tfj decay, but not at the correct 
kinematics for that decay. Likewise a signal is seen for El 
h c — > r] c decay in jpsigammaO method. We will discuss 
these transitions further elsewhere. 



From eq. (27) we can determine V(0) given a value 



for Vqq™ and a renormalisation factor, Z. For Z we 
use the fully nonperturbative vector form factor method 
described in Appendix |B 2| which normalises the local 
charm-charm vector current that we are using here by 
demanding that its form factor is 1 between identical 
mesons at q 2 = 0. This requires a non-staggered specta- 
tor quark and we use NRQCD for this. The determina- 
tion of Zff then needs the calculation of the form factor 
of the temporal component of the vector current between 
two £? c -like mesons (the mesons do not have to be real 
B c mesons) at rest. Zff can be determined with a sta- 
tistical uncertainty of 0.1% this way. Details are given in 
Appendix |B 2| 

The values for Z are given in Table [IX| of Appendix [B2 
and the values we use here are reproduced in Table VI 
along with our results for Vjjg™, V(0)/Z and the r) c and 
J/ip masses and energies. The table is divided into two 



with the upper results from the jpsigammaO method and 
the lower results from the etacgammaO method. The two 
methods give results for V(0)/Z in good agreement, but 
the jpsigammaO results are statistically more accurate. 
This is then our preferred method and the one that we 
will use for our final result. The agreement between the 
two methods to within the 2% statistical errors is a strong 
test of the control of discretisation errors in the HISQ 
formalism. 

Table |IX] also gives results that allow us to test to what 
extent V^(0) depends on m c and the precise tuning of q 2 to 
zero. On set 2 we have deliberately mistuned the c quark 
mass by 5% and see that it makes no significant difference 
to V(0) within our 2% statistical errors, q 2 is tuned to 
zero typically within our statistical errors of (lOMeV) 2 . 
On set 2 comparison between two different values of q 2 
shows no effect within our 1% statistical errors. We use 
the value closest to q 2 = in our fits below. These are 
both good tests of the robustness of our results to the 
tuning of parameters. 

Figure [7] shows our results for V^(0) plotted as a func- 
tion of the lattice spacing. To determine the physical 
value we use a fit similar to that for the hyperfinc split- 
ting and leptonic decay constant given in eq. [3j We sim- 
plify the fit slightly in dropping the tuning for the phys- 
ical c mass since our results in Table |VI| show negligible 
dependence on the c quark mass. We take the prior on 
the physical value to be 2.0(0.5) and allow for terms in 
(m c a) 21 up to i = 5. We take the prior on the leading 
(m c a) 2 term to be 0.0(3) since tree-level a 2 errors are re- 
moved in the HISQ action. We take linear and quadratic 
terms in 26xi + Sx s and allow a 2 dependence multiplying 
the linear term. 

The physical value for V(0) from the fit is 1.90(7) from 
the jpsigammaO method. The etacgammaO method gives 
a result in good agreement with a very similar error. The 
error is dominated by that from the extrapolation in the 
lattice spacing. In fact there is no visible lattice spacing 
dependence in our results and it could be argued that, 
in a transition from J /if) to rj c that probes relatively low 
momenta, the relevant scale for discretisation errors is 
well below m c . However, to be conservative, we allow 
discretisation errors to depend on (m c a) 2 and allow for 
multiple powers to appear. 

We have also tested extrapolations of V(0) to the phys- 
ical point using alternative definitions of the renormali- 
sation of the current. We get the same answer using Zft 
values taken from B c — > B c form factors with a heavier 
6 quark mass, as given in Appendix |B 2| We also get a 
result in good agreement if we use values for Z from Z cc 
given in Appendix |B 1| 

Our physical result for V^(0) is for a world that does 
not include electromagnetism, c-in-the-sea or allow for rj c 
annihilation. The effect of missing electromagnetism is 
similar to that for the decay constant and so we allow 
the same additional systematic error of 0.5%. We expect 
c-in the sea effects to be negligible, as for the decay con- 
stant. 7/ c annihilation affects the mass difference between 
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FIG. 8: A comparison of results for the charmonium hyper- 
fine splitting from lattice QCD with experiment. We show 
only results that include sea quarks and make use of multiple 
lattice spacings to derive a continuum value. The experimen- 
tal average [7] is given at the top, followed by the result for 
HISQ quarks from this paper. The Fermilab clover [29] and 
twisted mass [10] results follow. Neither of these lower two re- 
sults include an error for missing r) c annihilation effects. This 
error is the dominant error for our calculation. Here we show 
our error bar excluding this effect as a solid line and the total 
error including this effect as a dotted line. 



the J/ip and r\ c (as discussed in section III A I and there- 



fore affects the momentum of the r\ c that corresponds to 
q 2 = for this decay. Equivalently it means that the real 
q 2 = point corresponds to a non-zero q 2 in our calcu- 
lation. Since we allow an uncertainty in the r\ c mass of 
2.4 MeV (Table III) this corresponds to an uncertainty 
around q 2 — of 6 x 10 _6 GeV 2 , keeping the spatial mo- 
mentum fixed. From Table IVII we see that this would 
produce a negligible change in V(0), not visible beneath 
our statistical errors. In addition we can use informa- 
tion from [5] which used results at different q 2 values to 
extrapolate to q 2 — albeit in the quenched approxima- 
tion. The q 2 dependence gave a change in V(q 2 ) from 
V(0) of 20% when q 2 was lGeV 2 . From this it is clear 
that effects from a slight mistuning because of r\ c anni- 
hilation effects should be completely negligible. We take 
as our final result then: 



V(0) = 1.90(7)(1). 



(28) 



The complete error budget is given in Table III 



FIG. 9: A comparison of results for the decay constant of 
the J/ip from lattice QCD with experiment. We include only 
results that include sea quarks and make use of multiple lat- 
tice spacings to derive a continuum value. The experimental 
average [71 is given at the top, followed by the result for HISQ 
quarks from this paper. The twisted mass [ID] results follow. 



IV. DISCUSSION 

Figure [S] compares our result for the charmonium hy- 
perfine splitting to experiment and to that from other 
lattice QCD calculations. We only show results that have 
been obtained including sea quark effects and making use 
of multiple lattice spacing values to derive a physical con- 
tinuum result. Values are also given for different forms of 
the clover action in J30H33I but either at only one value 
of the lattice spacing or without giving a value from con- 
tinuum extrapolation. Some of these latter calculations 
obtain values well below experiment because of the large 
discretisation errors, particularly for the hyperfine inter- 
action, in the clover formalism. 

Our result agrees well with experiment and is more ac- 
curate than earlier values, especially since earlier values 
do not generally include any error for missing r\ c annihi- 
lation effects. 

Figure |9] similarly compares our result for fj/^, to that 
from twisted mass quarks including only u and d quarks 
in the sea [10] and to experiment (from eq. ([8])). Both 
lattice results agree well with experiment at the 2% level 
of accuracy achieved. Our value for fji^ gives a value 
for T{J/ip -> e+e-) of 5.48(16) keV using eq. (§]). 

Figure [10] shows the same comparison for the vector 
form factor at q 2 — 0, V(0), for J/ip — » 7^7 decay. Our 
result here using HISQ quarks and including u, d and s 
quarks in the sea agrees well, at the 4% level of accuracy 
achieved, with the result using twisted mass quarks and 
including only u and d sea quarks. 
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FIG. 10: A comparison of results for the vector form factor, 
V(0) for J/ip — > r/c^y from lattice QCD with experiment. We 
include only results that include sea quarks and make use of 
multiple lattice spacings to derive a continuum value. The 
experimental result [26] is given at the top, followed by the 
result for HISQ quarks from this paper. The twisted mass [TU] 
results follow. 



The value of V(0) extracted from the experimental 
branching fraction [55] is 1.7cr lower than the lattice num- 
bers where a is dominated by the 8% uncertainty from 
experiment. This situation is an improvement over that 
before the CLEO measurement [8j. However, it is clear 
that a more stringent test of QCD would be possible with 
a smaller experimental error for the J/ip — > r\ c ^ branch- 
ing fraction and this may become possible with BES III 
although it is a challenging mode [3U [35] . 

Our value for V(0) corresponds to a width for J/ip — > 
jr] c of 2.49(18)(7) keV using eq. (14). The first error is 



from our result and the second from the experimental er- 
ror in \q\ . Note that in using eq. ( 14 ) we put in the exper- 
imental masses for the J/ip and rj c . This is appropriate 
because these factors are kinematic ones and therefore 
should be taken to match the experiment. What we cal- 
culate in lattice QCD is V(0). In fact, as discussed above, 
we have good agreement between our results and experi- 



ment for M 



- M„ and so the kinematic factors would 



also be correct from lattice QCD. However, extra uncer- 
tainty would be introduced by using the lattice QCD re- 
sults and that is not necessary or appropriate. Our result 
for the decay width corresponds to a branching fraction 
for J/i/j T] c of 2.68(19)(11)%, where the first error is 
from our calculation and the second from experiment, 
including the experimental width of the J/ip- 

Figures [2j [3] and [7j which show our results as a function 
of lattice spacing, confirm that discretisation errors are 
small (although visible) for the HISQ formalism and that 



the approach to the continuum limit is well-controlled. 
This is discussed further in Appendix [C] where we com- 
pare the dependence on lattice spacing to that for twisted 
mass quarks [TU] , 



V. CONCLUSIONS 

We have given results for 3 key quantities associated 
with the J/ip meson from lattice QCD, for the first time 
including the effect of all three u, d and s quarks in the 
sea. The quantities are the mass difference with its pseu- 
doscalar partner, the r\ c meson, the decay constant and 
the vector form factor at q 2 = for J/ip — > rj c decay. 

Our first key result is for the J/ip decay constant. We 
obtain: 



f J/i> = 405(6) MeV, 



(29) 



leading to T(J/ip e+e") = 5.48(16) keV. This is to be 
compared to the experimental result of T(J/ip —> e + e~) 
= 5.55(14) keV [7J. We have therefore achieved a 4% test 
of lattice QCD from an electromagnetic decay rate (a 2% 
test from the decay constant), that does not suffer from 
CKM uncertainties. This is itself a stringent test of QCD 
and one for which lattice QCD is absolutely necessary; 
/ju could not be calculated this accurately with any 
other method. At the same time we are able to verify 
that the time- moments of the J/ip correlator agree as 
they should with results for the charm contribution to 
a(e + e~ — > hadrons) extracted from experiment. This is 
a test of QCD to better than 1.5%. 

Our fj/^ result is a critically important test for our 
calculations that determine the decay constants of the 
D s [35] and the D [35] [37] to a similar level of pre- 
cision. In particular it tests the HISQ formalism for c 
quarks |11| even more stringently than in the D and D s 
cases because the J/ip contains two c quarks and is a 
smaller meson, more sensitive to discretisation effects on 
the lattice. Combined with our earlier work on using the 
HISQ formalism for light quarks in f n and fx [T3I [351 158] . 
our result for fj/^ provides compelling evidence that we 
have the systematic errors in fo 3 and fo under control. 

We can improve our result for fj/^p further in future 
by using the vector form factor method of renormalisa- 
tion rather than the current-current correlator method. 
This will only be useful if improved experimental results 
become available. This is expected from BESIII 35 . 

A further test of QCD/Lattice QCD comes from the 
J/ip mass. We find: 



Mj/^ - M Vc = 116.5 ± 3.2 MeV 



(30) 



giving Mj/^p = 3.0975(32) (11) GeV where the second er- 
ror comes from the experimental average for M I)c [7j. 
Experiment gives Mj/^ = 3.0969 GeV. This is another 
strong test of lattice QCD, and indeed QCD, against ex- 
periment to be compared to that of the determination 
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of Md b [5] and Mjj [3] ■ The hyperfine splitting is a rel- 
atively small relativistic correction in the broader con- 
text of charmonium meson masses and the fact that we 
can do this well (with no free parameters) is because the 
HISQ formalism is such a highly improved relativistic 
formalism. This is underlined by a study of the meson 
dispersion relation (and associated 'speed of light') in 
Appendix [Cl In fact our error on Mj/^ is dominated by 



uncertainties from the effect of annihilation of the r\ c me- 
son to gluons, and it is important to pin these down more 
accurately. 

Our third result for the J/ip is that for its Ml radiative 
decay mode to the rj c . We find: 



r(J/V> jrjc) = 2.49 ± 0.19keV 



(31) 



to be compared to the current experimental value of 
1.84(30) keV [26]. The agreement is reasonably good, 
but the experimental error is large and the lattice QCD 
result would allow a much stronger test of QCD if this 
were reduced. This should be possible at BESIII [3"4"] . 
Since the error in our lattice QCD result is dominated by 
the continuum extrapolation it will be improved in cal- 
culations on superfine and ultrafine lattices as we have 
done for the decay constant, and 2% errors should also 
be achievable here. Again, this is only possible in lattice 
QCD. 

The J/ip — > 7?y c decay rate is another test of QCD, 
along with our leptonic decay constant, that is free of 
CKM uncertainties. It provides a validation of semilcp- 
tonic decay rate calculations for D and D s mesons [3"51 - 
|4"3"] . that also use HISQ quarks as well as a test of our 
techniques for nonperturbative current renormalisation 
that we are using for a range of semileptonic and radia- 
tive decays 
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FIG. 11: Difference in mass in MeV for different meson 
'tastes' for the r\ c and J/V> used here, plotted against the 
square of the lattice spacing. Red open circles show the mass 
difference between the local non-Goldstone and goldstone r/ c 
mesons. For the vector we have the mass difference between 
the local J/i/> meson and the 1-link ji CS> 'jijj vector (blue 
crosses) and the 1-link 707^ ® 7o7i7j vector (blue open tri- 
angles). The local J/tp meson is the lighter in both cases. 
Results are for sets 1, 2 and 4 from Tables [II] and |VT] Errors 
are statistical. 



Appendix A: Taste effects in staggered mesons 

Each staggered meson comes in 16 different tastes, 
most easily seen in terms of naive quark operators made 
with different point splittings: 



(Al) 



Here s is a 4-dimcnsional vector with or 1 in each com- 

(s) 

ponent. The different J„ operators are orthogonal to 
each other. To work out the corresponding staggered 
quark correlators we need the staggering matrix, Q (x). 
In our convention this is: 



n(*) = li(7/.) a 



(A2) 



with 74 = 7o- Then the connection between naive quark 
propagators, which carry a spin index, and staggered 
quark propagators, which do not, is: 

S F (x, y) = (V(a#(l/)>«. = 9(x, yMx)n\y). (A3) 

To work out the phases that appear in the correlator 
of a particular taste we then simply have to calculate 
spin traces over products of 51 and 7„ factors, see, for 
example, [llj . 

Here we use two different tastes for the r] c and three 
for the J/ip and they all have slightly different masses. 
The mass differences between the different tastes of a 
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jo) (i) -rzy 



J3T 



Set m c a Z(4) = Z cc Z(6) 



n=4 1.0 0.235 0.354 -0.187 0.0(5) 
n=6 1.0 0.246 0.460 0.198 0.0(5) 
n=8 1.0 0.253 0.563 0.511 0.0(5) 



TABLE VII: The perturbative series [45H49] for the ratio 



2 /r„ for different moments, n, in continuum QCD per- 



turbation theory, c^ 1 ' is the coefficient of ot-^ s , 

1.0 for all cases, by definition. for i > 3 were included in 
the determination of Z, allowing for a coefficient of 0.0 ± 0.5. 



„(o) 



given meson, however, vanish as a 2 a 2 - For the r\ c we 
use the Goldstone meson (in spin-taste notation this is 
the 75 ® 75 meson) and the local non-Goldstone meson 
(the 7075 <8> 7o75 meson), which are the first two mesons 
on the ladder of pseudoscalar tastes. We show the mass 
difference between them in Figure [IT] for coarse and fine 
lattices. The mass difference amounts on the coarse lat- 
tices to a little less than 10 MeV (for a 3 GeV particle) 
and clearly falls with a 2 as expected. For the J/ip we 
use the local vector (7^ <E) and two 1-link operators 
which have a point-splitting in an orthogonal direction 
to the polarisation (7^ (g) 7^- and 707^ ® 7i7j)- Note that 
these are not taste-singlet vectors. The mass difference 
between tastes for mesons of other J PC is typically much 
smaller than for pseudoscalars and that is clear here. Fig- 
ure [Tl] shows that the mass difference for the vectors is 
1-2 MeV on the coarse lattices and not resolvable on the 
fine lattices. 



In Section HI C| we showed results for J/ip — > 7?y c using 
different tastes of J/tp and rj c at the two ends of the 3- 
point function. No difference was seen in the vector form 
factor at q 2 = in the two cases, either on the coarse or 
fine lattices (within our statistical errors of 2%). This is 
another demonstration that taste effects are very small 
with HISQ quarks. 



Appendix B: Determining nonperturbative Z factors 
for local vector currents 



1. The current-current renormalization method 

Time-moments of lattice QCD correlators for zero- 
momentum heavyonium mesons can be compared very 
accurately |SD] to continuum QCD perturbation the- 
ory [4"5T44"5] developed for the analysis of the e + e~ an- 
nihilation cross-section. This has been used with pseu- 
doscalar meson correlators made with HISQ quarks to ex- 
tract c and b masses and a s to better than 1% |51| . These 
results used the goldstone pseudoscalar correlator, which 
is absolutely normalised because of the HISQ PCAC re- 
lation. Here we apply the same techniques to vector me- 
son correlators but use it to determine the renormalisa- 
tion factor, Z, required for the lattice vector current to 
match the continuum current. 

The time moments of our lattice QCD correlators are 



0.622 0.979(12) 
0.63 0.979(12) 
0.66 0.974(12) 
0.413 0.983(12) 
0.273 0.986(12) 
0.193 0.990(12) 



0.945(14) 
0.945(14) 
0.941(14) 
0.953(14) 
0.970(14) 
0.982(14) 



0.927(17) 
0.926(17) 
0.921(17) 
0.953(17) 
0.975(18) 
0.986(18) 



TABLE VIII: Renormalisation constants determined from 
the current-current correlator method on each configuration 
set used for the determination of fj/^,. The Z value we use 
is that from moment 4. The errors include an estimate of 
effects from a gluon condensate contribution, and unknown 
fourth order and higher terms in continuum perturbation the- 
ory. The errors are highly correlated between configuration 
sets (to better than 1% of the error). For set 2 we include 
both the tuned value of am c (0.63) and the heavier, detuned, 
value (0.66). Very little difference is seen between them. 



defined as: 



and 



(Bl) 



(B2) 



where t is a symmetrised version of t around the centre 
of the lattice, i.e. going forward in time, t runs from 
to L t /2 and then from —L t /2 + 1 to —1. The extra 
factor of (am c ) 2 in the pseudoscalar case is to make a 
correlator moment which is finite as a — > 0. For both 
correlators we expect the small n moments to behave 
perturbatively, since they probe small times. Then our 
match to continuum perturbation theory is: 



/■■ _ ^(W^.mK) 



cr = 



{am c (p)) n 



0((am c ) m ). 



(B3) 



where g n is the continuum QCD perturbation theory se- 
ries in the MS scheme [45 49J. For the vector correlator 
a Z factor is needed to multiply the lattice current and 
so: 



1 9n ("Mstrf^K) 
Z 2 {am c {fi)) n - 2 



+ 0((am c ) m ). (B4) 



Z is then a function of the bare lattice strong coupling 
constant at each lattice spacing and so this match must 
be performed separately on each ensemble. By taking the 
ratio of pseudoscalar and vector moments we can cancel 
the factors of the quark mass. In fact we also divide each 
moment by its tree-level value (calculated with the gluon 
fields set to 1) to reduce discretisation errors, i.e. instead 
of C n we use R n where 



Rn = 



C 



(0) 



(B5) 
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and also take the same ratio, calling it r n , in the contin- 
uum perturbation theory. Finally Z is given by: 



lK + 2/R V n 



' n+2 



I ' r, 



(B6) 



Table VII| gives the perturbative coefficients for the 



series in ctjjg(m) for r^ +2 /r^. This is known to 4- 
loops (ctg) and we include the possibility of unknown 
higher order terms (to 20th order) with prior values for 
the coefficients of ± 0.5. We take fx = m but in- 
cluding the possibility of higher order corrections means 
that the results are almost completely insensitive to /i. 
We also allow for gluon condensate contributions taking 
a s G 2 /n — ± 0.012 GeV 4 . This increases the errors on 
the determination of the Z values as n increases. 

Table |VIII| gives the results for Z on each ensemble 
and for moments 4, 6 and 8. The differences between the 
Z values on a given ensemble arise from discretisation 
errors. We take our final result for use in section IIIIBI 
from moment 4 (and so these numbers are reproduced in 
Table [n| since, even though discretisation errors fall as 
the moment number increases, the errors from the gluon 
condensate rise more steeply. The error in Z{4) is around 
1%. It is dominated by the uncertainty in higher orders 
in perturbation theory and so strongly correlated from 
one lattice spacing to the next. We have checked that we 
obtain the same final result for /w, using Z(6) (but with 
an error of 1.6% rather than 1.4%). 



2. The vector form factor method 

Vector currents can be normalised completely nonper- 
turbatively by requiring that the vector form factor at 
q 2 = between two identical mesons be 1 since this 
would be true for a conserved vector current. We use 
this to normalise the staggered taste-singlet 1-link vector 
current between two mesons made of staggered quarks 
in [U]. Here, however, we want to normalise the taste- 
nonsinglet local staggered vector current, and we cannot 
do this using a 3-point function made entirely of stag- 
gered quarks. We have to include a non-staggered (and 
non-doubling) spectator quark in order to remove the 
requirement for overall taste-cancellation in the 3-point 
function [52]. The Fermilab Lattice/MILC collaboration 
use this method 53 with a clover spectator quark to nor- 
malise a light staggered vector current for the nonpertur- 
bative part of their mixed perturbative/nonperturbative 
approach to normalising the clover-staggered operators 
for heavy-light meson decay constants. Here we use an 
NRQCD spectator quark to normalise the local charm 
staggered vector current. In fact we can use any NRQCD- 
like action for the spectator since it does not need to cor- 
respond to a physical quark. However, for simplicity we 
do use the NRQCD action developed for our NRQCD- 
light spectrum calculations [3|. 

To combine an NRQCD (or other non-staggered) quark 
with a staggered quark we convert it to a naive quark |52j , 




FIG. 12: The ratio of the 3-point function for the '£> c ' to 
'_B C ' vector form factor to the l B c ' 2-point function against the 
source-current separation, t/a. The top plot shows results for 
iriiO = 2.0 on fine lattices, set 4, for various values of T. Lines 
join the points for clarity. The shaded red band gives our fit 
result for the ground state matrix element. The lower plot 
shows results, also on set 4, for two different values of mta. 
The shaded red and blue bands show the fit results for the 
ground state matrix elements. 



re- instating the spin degree of freedom as in eq. ( A3 ) . A 



pseudoscalar meson correlator which, with two idential 
quarks, would simply be the sum over spins and colors of 
the squared modulus of the propagator, becomes in this 
case: 

CqS) ^J2 r Tr{G(x,y)n(y)g\x,y)n\x)} . (B7) 

The trace is over spin and color but separates since the 
staggered propagator, g(x, y), has no spin and Q no color. 
Then: 

CqS) = J2^ c {Tr s [^(x)G(x,y)n(y)}g^x,y)} . 

. (B8) 

This makes clear that we can transfer the fl matrices to 
the non-staggered quark propagator, G(x, y), and this al- 
lows us, as shown in |14j . to combine staggered and non- 
staggered quarks from random wall sources. We simply 
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make the source of the non-staggered quark the product 
of f2(y) (where y runs over a time slice at the source) with 
the random number at each y value in the random wall, 
convoluted with smearing functions if required. At the 
sink we multiply by W before tracing over spins to com- 
bine with the spinless staggered quark propagator. Here 
we develop this method further for 3-point functions. 

For the matrix element of the temporal charm-charm 
vector current between identical charmed pseudoscalar 
mesons, at rest, we have: 



(PcWt\P c 



2m P J+(0) 



(B9) 



and we can demand that /+(0) = 1, i.e. we can multiply 
the left-hand-side by Z to make this true. The temporal 
component of the vector current is the easiest to use for 
this purpose although the spatial component of the vector 
current is the one that we use for J/ip — » r/ c decay. For 
a relativistic action such as HISQ the renormalisation of 
the spatial and temporal components will be the same 
up to discretisation errors. 

For the 3-point function needed to evaluate the matrix 
element above we use an NRQCD heavy quark propagat- 
ing from to T and HISQ charm quarks from to t and 
T to t (see Fig. [5]) . The local staggered temporal vector 
current (70 ® 70 in spin-taste notation) is inserted at t 
and the restriction for staggered 3-point functions of an 
overall tasteless correlator is avoided by the spin content 
of the NRQCD propagator. 

The 3-point function is given by: 

C 3pt (0,t,T) = {-l) XT+VT+ZT+tT {-l) u (BIO) 

x Tr c {g(t, T)Tr s [ 7o Ot (T)G(T, O)fi(O)]^ (t, 0) } 

The g are HISQ c propagators and G is an NRQCD heavy 
quark propagator. The 70 factor comes from the local 
temporal vector current. The 3-point function is there- 
fore calculated in a very similar way to the 2-point func- 
tion. Q(x) multiplied by the random wall at the source 
time slice, 0, is used as the source for the NRQCD prop- 
agator. At time T this is multiplied by and 70, source 
and sink spins are set equal and summed over. This 
is then the source for the HISQ propagator from T to 
t which is finally combined with the HISQ propagator 
generated from the random wall at time slice 0. 

We calculate the 2-point and 3-point functions de- 
scribed above for several different NRQCD masses, m^a, 
and c quark masses, m c a, on the configuration sets 1, 2 
and 4. We also use several different values of T so that 
our fit benefits from both T and t dependence to improve 
the extraction of the ground-state masses, amplitudes 
and matrix elements. The 3-point and 2-point correla- 
tors arc fit simultaneously to the forms given in eqs ([I]) 
and ( 25 ) . We use the same priors as in section HI C| Note 
that for the 3-point fit we can now impose symmetry un- 
der interchange of the mesons at and T since they are 
the same. This means that the amplitudes V nn and V°° 
are square symmetric matrices and V no = V on . 



The key quantity that we extract from the fit is the 
ground state matrix element Vqq 1 . This is proportional 
to the vector matrix clement on the left-hand-side of 
eq. ( B9 1 . We can work out the constant of proportion- 



ality by matching our fit equations, eq (25) and ([I]) to 
the form expected by inserting a complete set of states 
in a continuum 3-point function. In this case factors of 
the mass of the meson 4 cancel and we find Vqq 1 = /+ (0) . 
Then the renormalisation factor we need is given by: 



Z 



V 00 



(Bll) 



The results for Z are taken from fits with 5 normal 
and 4 oscillating exponentials and given in Table |IX| We 
obtain very precise results for Z, with errors of 0.1%, 
without even using the full statistics available for each 
ensemble. The values are much more precise, for exam- 
ple, than for the implementation given in [53| (although 
their values are for light quarks rather than charm). 



Figure [12| shows the quality of our results through plots 
of the ratio of the 3-point to the 2-point function. Note 
that we fit the 3-point and 2-point functions simultane- 
ously and not just the ratio. It is convenient to plot the 
ratio, however, because the ground state contribution to 
this is simply Vqq 1 . Our fit allows us to include the ef- 
fect of excited states, both normal states and oscillating 
states. The presence of oscillating states is evident in the 



plots. The upper plot of Figure 12 compares results at 
different T values. All are included in the simultaneous 
fit. The lower plot shows results at different values of 
rnhd for a given m c a. 

The vector form factor method for determining Z is 
completely nonperturbative. It will therefore be subject 
to errors coming from lattice QCD in the form of dis- 
cretisation errors. These mean, for example, that the 
Z factor at a particular value of the lattice spacing is 
not completely independent of the mass of the NRQCD 
spectator quark. From our results in Table |IX| we see 
that changing m^a from 2.8 to 1.5 on coarse set 2 (cor- 
responding to change of almost a factor of two) causes a 
2% change in Z. On fine set 4 the sensitivity is reduced 
to a change of 0.2%, not significant within our statistical 
errors, when mha is changed from 2.0 to 1.5, a change 
in mass of 30%. Since the change in lattice spacing be- 
tween the coarse and fine sets is a factor of 1.4, pairs of 
arrih values on coarse and fine lattices that correspond 
to approximately the same physical mass are (2.8, 2.0) 
and (2.0, 1.5). We will take our central result for Zff 
from the Z values corresponding to using (2.0, 1.5) and 
use (2.8, 2.0) to check systematic errors coming from Z. 



We see in section III C that we get the same answer from 
both sets of Z values. 



4 For a meson containing an NRQCD quark the energy obtained 
from the 2-point and 3-point fits is not its mass. However that 
is irrelevant here since it cancels. 
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Set TVcfg x N t 


T arrih 


am c 


e aE Pc VqT Z = Z ff 


1 


450 x 4 


20,21 2.0 


0.622 


-0.221 0.9630(2) 1.0104(12) 0.9896(11) 


2 


408 x 4 


20,21,24 2.8 
20,21,24 2.0 
20,21 1.5 
20,21 2.0 


0.63 
0.63 
0.63 
0.66 


-0.226 1.0239(2) 1.0220(10) 0.9784(9) 
-0.226 0.9719(2) 1.0106(8) 0.9894(8) 
-0.226 0.9311(3) 1.0026(14) 0.9974(14) 
-0.244 0.9994(3) 1.0138(18) 0.9863(17) 


4 


322 x 4 


24,25,30 2.0 
24,25,30 1.5 


0.413 
0.413 


-0.107 0.6454(2) 0.9966(14) 1.0033(14) 
-0.107 0.5939(2) 0.9950(10) 1.0049(10) 



TABLE IX: The Z factors (column 9) obtained from the vector form factor method on different configurations sets (column 
1) and for different NRQCD quark masses, mha (column 4), and c quark masses, m c a (column 5) with e factor (column 6). 
The m c a values are those used in our calculation of J/ip — > r\ c ^ described in section |ill C| the m^a values are arbitrary since 
they correspond to the spectator quark. Z is given by the inverse of the fit parameter Voo™ given in column 8. Column 2 gives 
the number of configurations used in the calculation and the number of time sources for the origin, 0. The T values used are 
given in column 3. Column 7 gives the energy of the NRQCD-c meson obtained from the fit. This is not equal to the mass 
because there is an energy offset in NRQCD. 



Set m c a aM n 



\ap\ 



aE v 



2 0.63 1.80851(4) 0.52880 1.88286(12) 0.9814(15) 

0.35000 1.84123(5) 0.9748(5) 
0.20000 1.81923(4) 0.9715(6) 

4 0.413 1.28042(4) 0.37486 1.33352(6) 0.9878(8) 
0.20000 1.29575(4) 0.9873(7) 

TABLE X: Rest masses («M,J and energies (aE Vc ) at non- 
zero momentum \ap\ for the Goldstone r\ c meson on sets 2 
(coarse) and 4 (fine). The rest masses differ slightly from 
those in Table [IT] because they come from independent fits; 
on set 4 we have higher statistics here. The zero and non- 
zero momentum correlator are fitted s imu ltaneously and the 
speed of light, c 2 , extracted using eq. (CI I. 



Table |IX| also shows results for the deliberately mis- 
tuned c quark mass of 0.66 on set 2 with which we test 
the c quark mass dependence of our J /%p — > r\ c form fac- 
tor. A barely significant change in Z is seen between this 
value and that for the tuned mass of 0.63. We also see 
no significant difference in Z as we change the sea light 
quark masses between set 1 and set 2. 

Discretisation errors also mean that our results for Zff 
here do not have to agree exactly with our earlier re- 
sults for Z cc . Z cc has a much larger error coming from 
unknown higher order terms in continuum perturbation 
theory. As a result of this, Zff and Z cc do agree at the 
level of la. They also agree within errors with the lattice 
QCD perturbation theory for this renormalisation, which 
has a small negative contribution at 0(a s ) |54) . 

The Z values calculated here form part of a programme 
to normalise nonperturbatively a range of staggered cur- 
rents for a variety of weak semileptonic and electromag- 
netic radiative decay rates [32] . 



Appendix C: Discretisation errors 

Finally we discuss discretisation errors. For relativistic 
formalisms the scale of discretisation errors can be set 
by the quark mass when this is larger than Aqcd and 
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coarse m c a=0.63 i — e- 
fine m a^).413 i — x— 
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(pa) 2 
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FIG. 13: The 'speed-of-light', c , calculated from zero and 
finite momentum r\ c correlators using HISQ c quarks. The 
top figure shows c 2 as a function of the square of the spatial 
momentum for coarse lattices, set 2, where m c a = 0.63 (grey 
open circles) and for fine lattices, set 4, where m c a — 0.413 
(pink crosses). The lower figure shows the resulting values of 
c 2 as ffa 2 — > as a function of (m c a) 2 . The dashed straight 
line is drawn to guide the eye. 
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FIG. 14: The rest (static) and kinetic masses for the r\ c me- 
son compared between the HISQ formalism (this paper) and 
the twisted mass formalism (Figure 3 from |10|). and plotted 
against the square of the lattice spacing. The rest mass is 
given by open circles and the kinetic mass by open triangles, 
red for twisted mass and blue for HISQ. Errors include the 
full lattice spacing error on each point. 



FIG. 16: fj/ij, in GeV is plotted against a 2 in fm 2 for the 
HISQ formalism (this paper, blue open circles) and for twisted 
mass ([10]. red open squares). For the twisted mass case the 
heaviest and lightest sea quark masses are plotted at each 
value of the lattice spacing. Only statistical errors are shown 
for the twisted mass results. For the HISQ results we show the 
raw data from Table ITT1 with statistical and uncorrelated lat- 
tice spacing errors. There is an additional error of 1.3% from 
correlated lattice spacing and Z factor uncertainties. The 
black cross is the experimental result from the average lep- 
tonic width [7], offset slightly from the origin for clarity. 
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FIG. 15: Rj/ii> — 1 plotted against a 2 in fm 2 for three 
different quark formalisms for the c quark: HISQ (this pa- 
per, blue open circles), twisted mass ([ID], red open squares) 
and Fermilab clover ( 29 , green open triangles, showing re- 
sults on the two finest lattices only). Rj/^ is Mj^/M Vc so 
Rj/i> ~ 1 — AMhy P /M Vc . For the twisted mass and Fermilab 
clover results the heaviest and lightest sea quark masses are 
plotted at each value of the lattice spacing. Only statistical 
errors are shown. Additional errors from (twice) the lattice 
spacing error amount to 2% for HISQ and Fermilab clover and 
4-7% for twisted mass. The black cross is the experimental 
average 7 1, offset slightly from the origin for clarity. 



1.2 1 ' ' 1 ' ' ' ' — 
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FIG. 17: The vector form factor for J/ij) — > rj c decay at 
q 2 = 0, V(0), is plotted against a 2 in fm 2 for the HISQ for- 
malism (this paper, jpsigammaO method, blue open circles) 
and for twisted mass ([ID], red open squares). For the twisted 
mass case the heaviest and lightest sea quark masses are plot- 
ted at each value of the lattice spacing. Errors include sta- 
tistical errors and uncertainties in the Z factors. The black 
cross is the experimental result from the rate for J/tp — > r/cj 
decay [26], offset slightly from the origin for clarity. 
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so they must be monitored closely when working with c 
quarks. 

One simple way to do this is through study of the 
energy of mesons at non-zero spatial momentum. Be- 
cause the HISQ formalism is a relativistic formalism we 
determine meson masses as the result of fitting zero- 
momentum meson correlators as described in section HT1 
For heavy quarks discretisation errors mean that this 
mass, known as the 'rest' or 'static' mass, can differ from 
the mass that controls the momentum-dependence of the 
energy at non-zero spatial momentum. This latter mass 
is known as the 'kinetic mass'. An equivalent statement 
is that the square of the speed-of-light, c 2 , differs from 1, 
where 111] 



In Figure 14 we compare the rest and kinetic masses for 



c 2 (p) 



E 2 (p) 



P 2 



(CI) 



For r/c mesons we are able to determine c 2 very accu- 
rately with HISQ quarks when we have 0(10,000) cor- 
relators as here. We fit zero and non-zero momentum 
simultaneously to the form given in eq. ([I]) (although the 
zero momentum correlators have no oscillating compo- 
nent). From the fit we obtain the ground state mass, 
M , in the zero momentum case and the ground state 
energy, Eq, in the non-zero momentum case. The simul- 
taneous fit allows us to take correlations into account to 
improve the error in c 2 . Results are given in Table |x[ c 2 
is within 3% of 1 but we can distinguish it from 1 and 
we can see that it depends on the spatial momentum. 



This is shown in the top plot of Figure 13 for the coarse 
lattices, set 2, and the fine lattices, set 4. All tree-level 
a errors are removed in HISQ but c 2 can depend lin- 
early on aPp 2 through 0(a s ) corrections. We see that 
the slope of c 2 with a 2 ]] 2 is much smaller on the fine lat- 
tices than the coarse, as m c a is reduced. This plot can 
be compared with earlier results in Figure 9 of |55j . al- 
though note that those results show a jump as the lattice 
momentum changes from (1,1,1) to (2,0,0). This results 
from a rotationally-noninvariant discretisation error not 
evident in our results because of the use of the phased 
boundary condition of eq. ( 16 1 to fix the momentum di- 



rection and simply change its magnitude. 

From the top plot of Figure [13] we can determine the 
value of c 2 at cPp 2 = 0. The results from the coarse and 
fine lattices are then shown in the lower plot to lie on 
a straight line as a function of (m c a) 2 . The small slope 
is compatible with the (m c a) 2 dependence also being an 
0(a s ) effect. The straight line clearly goes through 1 as 
m c a as it must. 

Another way to look at these discretisation errors is 
to compare the rest and kinetic masses for the r] c . The 
kinetic mass is given by [56 : 



Mhin. — 



\p\ 2 - (A£) 2 
2AE 



(C2) 



where AE = E(p) — M rest . We also have c 2 = 
Mrest/M^in [TT] . In the absence of errors, for a relativis- 



HISQ 77 c mesons using M kin = M rest /((P(p^ = 0)) de- 
termined from Figure [l3j The rest and kinetic masses 
differ by 3% on the coarse lattices and 1.3% on the fine 
lattices. 

Figure [14] also compares results from the twisted mass 
formalism, from [10] . For that formalism there is a signif- 
icantly larger difference between rest and kinetic masses 
for the f] c meson, amounting to 35% on the coarsest lat- 
tice spacing. The twisted mass formalism has tree-level 
a 2 errors, so a larger effect would be expected. The rest 
and kinetic masses agree in the continuum limit, as they 
must. 

Figure [15] compares the hyperfine splitting in charmo- 
nium in units of the r\ c mass as a function of the square 
of the lattice spacing for the HISQ and twisted mass for- 
malisms. The quantity plotted is Rj/^, — 1 where Rj/^p 
is defined in [TO] as Mj/^/M Vc and results are given on 
each of their ensembles. The HISQ results are taken from 
Table llll We also show results from the Fermilab clover 
method [53] on their two finest sets of ensembles which 
correspond to the coarse and fine lattices used here. 

All the results tend to the same continuum value which 
agrees with experiment. Much larger discretisation ef- 
fects are visible for the twisted mass formalism than for 
HISQ. These are compatible with the tree level a 2 errors 
expected in that formalism, and with a mass scale of ap- 
proximately 2 GeV (i.e. m c ). These tree level errors are 
removed in the HISQ formalism, but a s a 2 errors remain. 
These seem to be small for this quantity. The Fermilab 
clover discretisation errors, although in principle a s a, are 
also relatively small over this range of a. 

Discretisation effects in the hyperfine splitting can also 
enter through tuning of the c quark mass, because the hy- 
perfine splitting is very sensitive to this. As discussed 
earlier in this section, there can be significant differ- 
ences between rest and kinetic masses for mesons made of 
heavy quarks, and either can be used to tune the quark 
mass. Both relativistic formalisms, HISQ and twisted 
mass, use the rest mass. The Fermilab formalism uses 
the kinetic mass. 

For fj/$ the difference in discretisation effects between 
the HISQ and twisted mass formalisms is not as large. 



tic formalism, we should have M r 



M k 



(i.e. c 



1 



This is shown in Figure 16 Again, answers in agreement 
are obtained in the continuum limit. For the moments 
of the vector correlator our results for G\ (FigjZJ), which 
show very little dependence on the lattice spacing, can be 
compared to those for the twisted mass formalism in [ST] , 
where somewhat larger discretisation effects are visible. 

Finally in Figure [TjJ we compare results for the vector 
form factor for J/ip — > r\ c decay, V(0). Large discretisa- 
tion effects are evident in the twisted mass results. Once 
again discretisation effects in the HISQ results are small. 
Agreement in the continuum limit is again clear, however. 

The HISQ results shown here give a further and more 
testing demonstration than that of [TTJ [36] of how small 
the discretisation errors for the HISQ action are, even for 
a quark as heavy as charm. 
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